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We propose a local real-space formulation for orbital-free DFT with density dependent kinetic 
energy functionals and a unified variational framework for computing the conhgurational forces 
associated with geometry optimization of both internal atomic positions as well as the cell geometry. 
The proposed real-space formulation, which involves a reformulation of the extended interactions 
in electrostatic and kinetic energy functionals as local variational problems in auxiliary potential 
helds, also readily extends to all-electron orbital-free DFT calculations that are employed in warm 
dense matter calculations. We use the local real-space formulation in conjunction with higher-order 
hnite-element discretization to demonstrate the accuracy of orbital-free DFT and the proposed 
formalism for the Al-Mg materials system, where we obtain good agreement with Kohn-Sham DFT 
calculations on a wide range of properties and benchmark calculations. Finally, we investigate the 
cell-size effects in the electronic structure of point defects, in particular a mono-vacancy in Al. 
We unambiguously demonstrate that the cell-size effects observed from vacancy formation energies 
computed using periodic boundary conditions underestimate the extent of the electronic structure 
perturbations created by the defect. On the contrary, the bulk Dirichlet boundary conditions, 
accessible only through the proposed real-space formulation, which correspond to an isolated defect 
embedded in the bulk, show cell-size effects in the defect formation energy that are commensurate 
with the perturbations in the electronic structure. Our studies suggest that even for a simple defect 
like a vacancy in Al, we require cell-sizes of ~ 10^ atoms for convergence in the electronic structure. 


I. INTRODUCTION 

Electronic structure calculations have played an im¬ 
portant role in understanding the properties of a wide 
range of materials systemsi. In particular, the Kohn- 
Sham formalism of density functional theory^ii^ has been 
the workhorse of ground-state electronic structure calcu¬ 
lations. However, the Kohn-Sham approach requires the 
computation of single-electron wavefunctions to compute 
the kinetic energy of non-interacting electrons, whose 
computational complexity typically scales as 0{N^) for 
an A^-electron system, thus, limiting standard calcula¬ 
tions to materials systems containing few hundreds of 
atoms. While there has been progress in developing 
close to linear-scaling algorithms for the Kohn-Sham ap¬ 
proach^, these are still limited to a few thousands of 
atoms, especially for metallic systems^. The orbital-free 
approach to DFT—, on the other hand, models the ki¬ 
netic energy of non-interacting electrons as an explicit 
functional of the electron density, thus circumventing the 
computationally intensive step of computing the single¬ 
electron wavefunctions. Further, the computational com¬ 
plexity of orbital-free DFT scales linearly with the system 
size as the ground-state DFT problem reduces to a min¬ 
imization problem in a single field—the electron density. 
The past two decades has seen considerable progress in 
the development of accurate models for orbital-free ki¬ 
netic energy functionals^A^^ and, in particular, for sys¬ 
tems whose electronic-structure is close to a free electron 
gas (for e.g. Al, Mg). Also, orbital-free DFT calculations 
are being increasingly used in the simulations of warm 
dense matter where the electronic structure is close to 
that of a free electron gas at very high temperaturesi^^— . 
As the reduced computational complexity of orbital-free 


DFT enables consideration of larger computational do¬ 
mains, recent studies have also focused on studying ex¬ 
tended defects in Al and Mg, and have provided impor¬ 
tant insights into the energetics of these defects2ir.2£. 

The widely used numerical implementation of orbital- 
free DFT is based on a Fourier space formalism using a 
plane-wave discretization22i2^. A Fourier space formula¬ 
tion provides an efficient computation of the extended in¬ 
teractions arising in orbital-free DFT—electrostatics and 
kinetic energy functionals—through Fourier transforms. 
Further, the plane-wave basis is a complete basis and pro¬ 
vides variational convergence in ground-state energy with 
exponential convergence rates. However, the Fourier 
space formulations are restricted to periodic geometries 
and boundary conditions that are suitable for perfect 
bulk materials, but not for materials systems contain¬ 
ing extended defects. Also, the extended spatial nature 
of the plane-wave basis affects the parallel scalability of 
the numerical implementation and is also not suitable for 
multi-scale methods that rely on coarse-graining. In or¬ 
der to address the aforementioned limitations of Fourier 
space techniques, recent efforts have focussed on devel¬ 
oping real-space formulations for orbital-free DFT and 
numerical implementations based on finite-element^Sri^i 
and finite difference discretizations^^Ti^i. 

In the present work, we build on these prior efforts to 
develop an efficient real-space formulation of orbital-free 
DFT employing the widely used non-local Wang-Govind- 
Carter (WGG)ii kinetic energy functional. As in prior 
efforts22i^j we reformulate the extended interactions in 
electrostatics and the non-local terms in the WGC ki¬ 
netic energy functionals as local variational problems in 
auxiliary potential fields. However, the proposed refor¬ 
mulation of electrostatic interactions is notably differ- 
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ent from previous works, and enables the evaluation of 
variational configurational forces corresponding to both 
internal atomic relaxations as well as external cell relax¬ 
ation under a single framework. Further, the proposed 
formulation naturally extends to all-electron orbital-free 
DFT calculations of warm dense matter— di. In the pro¬ 
posed real-space formulation, the ground-state orbital- 
free DFT problem is reformulated as an equivalent sad¬ 
dle point problem of a local functional in electron den¬ 
sity, electrostatic potential and the auxiliary potential 
fields (kernel potentials) accounting for the extended in¬ 
teractions in the kinetic energy functional. We employ 
a higher-order finite-element basis to discretize the for¬ 
mulation, and demonstrate the optimal numerical con¬ 
vergence of both the ground-state energy and configura¬ 
tional forces with respect to the discretization. Further, 
we propose an efficient numerical approach to compute 
the saddle point problem in electron density, electrostatic 
potential and kernel potentials by expressing the saddle 
point problem as a fixed point iteration problem, and 
using a self-consistent field approach to solve the fixed 
point iteration problem. 

We subsequently investigate the accuracy and transfer- 
ability of the proposed real-space formulation of orbital- 
free DFT for A1 and Mg materials systems. To this end, 
we compute the bulk properties of Al, Mg and Al-Mg in- 
termetallics, and compare it with Kohn-Sham DFT. As 
orbital-free DFT only admits local pseudopotentials, the 
Kohn-Sham DFT calculations are conducted using both 
local and non-local psedupotentials. Our studies indi¬ 
cates that the bulk properties computed using orbital- 
free DFT for Al, Mg and Al-Mg intermetallics are in good 
agreement with Kohn-Sham DFT. We further investigate 
the accuracy of orbital-free DFT by computing the inter¬ 
atomic forces in Al and Mg, which are also in good agree¬ 
ment with Kohn-Sham DFT calculations. Our studies 
demonstrate that orbital-free DFT is accurate and trans¬ 
ferable across a wide range of properties for Al, Mg and 
Al-Mg intermetallics, and can be used to study properties 
of these materials systems that require computational do¬ 
mains that are not accessible using Kohn-Sham DFT. For 
instance, in the present study we computed the forma¬ 
tion energy of /?' Al-Mg alloy containing 879 atoms in a 
unit cell employing the proposed real-space formulation 
of orbital-free DFT, but the same system was found to 
be prohibitively expensive using Kohn-Sham DFT. 

We finally investigate the cell-size effects in the elec¬ 
tronic structure of point defects, in particular a mono¬ 
vacancy in Al. Prior studies using Fourier-based for¬ 
mulations of orbital-free DFT have suggested that the 
formation energy of a mono-vacancy in Al is well con¬ 
verged by 108-256 atom cell-sizes^^. However, coarse¬ 
grained real-space calculations have suggested that much 
larger cell-sizes of the order of 1,000 atoms are required 
for convergence of vacancy formation energies^, which 
was also supported by asymptotic estimates^. In or¬ 
der to understand the underpinnings of this discrepancy, 
we use the finite-element discretized real-space formula¬ 


tion of orbital-free DFT and compute the vacancy for¬ 
mation energy using two boundary conditions: (i) pe¬ 
riodic boundary conditions, equivalent to Fourier-space 
based formulations; (ii) bulk Dirichlet boundary condi¬ 
tions, where the perturbations in the electronic struc¬ 
ture arising due to the vacancy vanishes on the boundary 
of the computational domain. Our study suggests that 
while the vacancy formation energy is well converged by 
108 atom cell-size using periodic boundary conditions, 
the electronic fields are not well-converged by this cell- 
size. On the other hand the bulk Dirichlet boundary 
conditions show well converged formation energy as well 
as electronic fields by cell sizes of ~ 1,000 atoms, which is 
consistent with prior real-space calculations. This study 
reveals that while periodic boundary conditions show a 
superior convergence in formation energies due to the 
variational nature of the formalism, the true cell-size ef¬ 
fects which also measure convergence of electronic fields 
are provided by the bulk Dirichlet boundary conditions. 
We note that the proposed real-space formulation with 
finite-element discretization are crucial to employing bulk 
Dirichlet boundary conditions, which enable the study of 
isolated defects in bulk. 

The remainder of the paper is organized as follows. 
Section II provides a description of the orbital-free DFT 
problem. Section III presents the proposed real-space 
formulation of the orbital-free DFT problem, the con¬ 
figurational forces associated with structural relaxations, 
and the finite-element discretization of the formulation. 
Section IV discusses the numerical implementation of the 
formulation and presents an efficient numerical approach 
for the solution of the saddle point real-space variational 
problem. Section V presents the numerical convergence 
results of the finite-element discretization of the real- 
space formulation, the accuracy and transferability of the 
real-space orbital-free DFT formalism for Al-Mg materi¬ 
als system, and the study of the role of boundary con¬ 
ditions on the cell-size effects in electronic structure cal¬ 
culations of point defects. We finally conclude with a 
summary and outlook in Section VI. 


II. ORBITAL-FREE DENSITY FUNCTIONAL 
THEORY 

The ground-state energy of a charge neutral materials 
system containing M nuclei and N valence electrons in 
density functional theory is given byiii 

R-) = Ts{p) + Exc{p) + Eh{p)+ Eextip, R-)+7?zz(R) , 

( 1 ) 

where p denotes the electron-density and R = 
{Ri, R2,..., Rm} denotes the vector containing the po¬ 
sitions of M nuclei. In the above, Tg denotes the kinetic 
energy of non-interacting electrons, E^c is the exchange- 
correlation energy, Eh is the Hartree energy or classical 
electrostatic interaction energy between electrons, E^xt 
is the classical electrostatic interaction energy between 
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electrons and nuclei, and E^z denotes the electrostatic 
repulsion energy between nuclei. We now discuss the var¬ 
ious contributions to the ground-state energy, beginning 
with the exchange-correlation energy. 

The exchange-correlation energy, denoted by E^c, in¬ 
corporates all the quantum-mechanical interactions in 
the ground-state energy of a materials system. While 
the existence of a universal exchange-correlation energy 
as a functional of electron-density has been established by 
Hohenberg, Kohn and Sharo^i^, its exact functional form 
has been elusive to date, and various models have been 
proposed over the past decades. For solid-state calcula¬ 
tions, the local density approximation (LDA) ^^i^^ and 
the generalized gradient approximatio n^^! have been 
widely adopted across a range of materials systems. In 
particular, the LDA exchange-correlation energy, which 
is adopted in the present work, has the following func¬ 
tional form: 


£'xc{p) 


= / e> 


,(p)p(x)dx 


( 2 ) 


where e^^{p) = ea;(p) -I- £c(p), and 

= (3) 


e^{p) = J (l-l-/3li/(r-„)-1-/32 r,) 

[ A\ogrs+B + Crs\ogrs + Drs < 1 , 

( 4 ) 

and Tg = (3/47rp)^/^. In the present work, we use the 
Ceperley and Alder constants^ in equation ()^. 

The last three terms in equation o represent elec¬ 
trostatic interactions between electrons and nuclei. The 
Hartree energy, or the electrostatic interaction energy be¬ 
tween electrons, is given by 


Eh{p) = 2 


P(x)p(xO 

|x-x'| 


dx dx'. 


( 5 ) 


The interaction energy between electrons and nuclei, in 
the case of local pseudopotentials that are adopted in the 
present work, is given by 


Eext^P: H-) — p{p^')^ext{p^: H-) dx 

= J2f P(x)14pi(|x-R/|)dx, (6) 

j '' 

where denotes the pseudopotential corresponding to 
the nucleus, which, beyond a core radius is the 
Coulomb potential corresponding to the effective nuclear 
charge on the nucleus. The nuclear repulsive energy 
is given by 


Ezz{^) 


I j.j^i 


ZiZj 

IR/ - R,/| ’ 


( 7 ) 


where Zj denotes the effective nuclear charge on the 
nucleus. The above expression assumes that the core ra¬ 
dius of the pseudopotential is smaller than internuclear 
distances, which is often the case in most solid-state ma¬ 
terials systems. We note that in a non-periodic setting, 
representing a finite atomic system, all the integrals in 
equations are over and the summations in 

equations (ED-© include all the atoms. In the case of 
an infinite periodic crystal, all the integrals over x in 
equations ©-© are over the unit cell whereas the inte¬ 
grals over x' are over Similarly, in equations ®-o, 
the summation over / is on the atoms in the unit cell, 
and the summation over J extends over all lattice sites. 
Henceforth, we will adopt these notions for the domain 
of integration and summation. 

The remainder of the contribution to the ground-state 
energy is the kinetic energy of non-interacting electrons, 
denoted by Tg, which is computed exactly in the Kohn- 
Sham formalism by computing the single-electron wave- 
functions (eigenfunctions) in the mean-field^. The con¬ 
ventional solution of the Kohn-Sham eigenvalue problem, 
which entails the computation of the lowest N eigenfunc¬ 
tions and eigenvalues of the Kohn-Sham Hamiltonian, 
scales as 0{N^) that becomes prohibitively expensive 
for materials systems containing a few thousand atoms. 
While efforts have been focused towards reducing the 
computational complexity of the Kohn-Sham eigenvalue 
problem^, this remains a significant challenge especially 
in the case of metallic systems. In order to avoid the 
computational complexity of solving for the wavefunc- 
tions to compute Tg, the orbital-free approach to DFT 
models the kinetic energy of non-interacting electrons as 
an explicit functional of electron density^. These models 
are based on theoretically known properties of Tg for a 
uniform electron gas, perturbations of uniform electron 
gas, and the linear response of uniform electron gas^^— . 
As the orbital-free models for the kinetic energy func¬ 
tional are based on properties of uniform electron gas, 
their validity is often limited to materials systems whose 
electronic structure is close to a free electron gas, in par¬ 
ticular, the alkali and alkali earth metals. Further, as 
the orbital-free approach describes the ground-state en¬ 
ergy as an explicit functional of electron-density, it limits 
the pseudopotentials calculations to local pseudopoten¬ 
tials. While these restrictions constrain the applicability 
of the orbital-free approach, numerical investigation s ^ ^4° 
indicate that recently developed orbital-free kinetic en¬ 
ergy functionals and local pseudopotentials can provide 
good accuracy for A1 and Mg, which comprise of tech¬ 
nologically important materials systems. Further, there 
are ongoing efforts in developing orbital-free kinetic en¬ 
ergy models for covalently bonded systems and transition 
metals^i^. 

In the present work, we restrict our focus to the Wang- 
Goving-Carter (WGC) density-dependent orbital-free ki¬ 
netic energy functionalii, which is a widely used kinetic 
energy functional for ground-state calculations of mate¬ 
rials systems with an electronic structure close to a free 
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electron gas. In particular, the functional form of the 
WGC orbital-free kinetic energy functional is given by 

Ts{p)=Cf J p^/^(x)dx-f i J |VA/p(x)pdx-f Tk(p) 

( 8 ) 

where 

Tk{p)=Cf j j p“(x) i4:(,^^(x, x'), |x - x'I) p^(x') dxdx', 
^^(x,x ) = - - -j , kF{x) = (37r>(x)) . 


In equation ([ 8 ]), the first term denotes the Thomas-Fermi 
energy with Cf = and the second term de¬ 

notes the von-Weizsdcker correction^. The last term de¬ 
notes the density dependent kernel energy, , where the 
kernel K is chosen such that the linear response of a uni¬ 
form electron gas is given by the Lindhard response^. 
In the WGC functionalii, the parameters are chosen to 
be {a, 13} = {5/6 -f v^/ 6 , 5/6 - ^5/6} and 7 = 2.7. For 
materials systems whose electronic structure is close to 
a free-electron gas, the Taylor expansion of the density 
dependent kernel about a reference electron density (po)j 
often considered to be the average electron density of the 
bulk crystal, is employed and is given by 


iF(f^(x,x'), |x 


x'l) =iFo(|x-x'|) -HiFi(|x-x'|)(Ap(x) -f Ap(x')) -f ^ATndx - x'|)((Ap(x))2-p (Ap(x'))^) 
-I-A:i 2 (|x - x'|)Ap(x)Ap(x')-f ... . 


In the above equation, Ap(x) = p(x) — po a-nd the density 
independent kernels resulting from the Taylor expansion 
are given by 


ATodx-x'l) 
ATidx-x'l) = 
ATiidx-x'l) = 
ATiadx-x'l) = 


|x- xd) 

9p(x) 

dp‘^{x) 

^^^(^7,|X-Xd) 

dp{x)dp{x') 


P=Po 


P=Po 


P=PO 


P=PO 


( 10 ) 


Numerical investigations have suggested that the Taylor 
expansion to second order provides a good approxima¬ 
tion of the density dependent kernel for materials systems 
with electronic structure close to a free electron ga o^^'^^ . 
In particular, in the second order Taylor expansion, the 
contribution from K 12 has been found to dominate con¬ 
tributions from Kii. Thus, in practical implementations, 
often, only contributions from K 12 in the second order 
terms are retained for computational efficiency. 


III. REAL-SPACE FORMULATION OF 
ORBITAL-FREE DFT 


In this section, we present the local variational real- 
space reformulation of orbital-free DFT, the configura¬ 
tional forces associated with internal ionic relaxations 
and cell relaxation, and the finite-element discretization 
of the formulation. 


A. Local real-space formulation 

We recall that the various components of the ground- 
state energy of a materials system (cf. section |n| are 
local in real-space, except the electrostatic interaction 
energy and the kernel energy component of the WGG 
orbital-free kinetic energy functional that are extended 
in real-space. Gonventionally, these extended interac¬ 
tions are computed in Fourier space to take advantage 
of the efficient evaluation of convolution integrals using 
Fourier transforms. For this reason, Fourier space for¬ 
mulations have been the most popular and widely used 
in orbital-free DFT calculations^iiS^. However, Fourier 
space formulations employing the plane-wave basis re¬ 
sult in some significant limitations. Foremost of these is 
the severe restriction of periodic geometries and bound¬ 
ary conditions. While this is not a limitation in the study 
of bulk properties of materials, this is a significant limi¬ 
tation in the study of defects in materials. For instance, 
the geometry of a single isolated dislocation in bulk is 
not compatible with periodic geometries, and, thus, prior 
electronic structure studies have mostly been limited to 
artificial dipole and quadrapole arrangements of dislo¬ 
cations. Further, numerical implementations of Fourier- 
space formulations also suffer from limited scalability on 
parallel computing platforms. Moreover, the plane-wave 
discretization employed in a Fourier space formulation 
provides a uniform spatial resolution, which is not suit¬ 
able for the development of coarse-graining techniques— 
such as the quasi-continuum method^^— that rely on an 
adaptive spatial resolution of the basis. 

We now propose a real-space formulation that is de¬ 
void of the aforementioned limitations of a Fourier space 
formulation. The proposed approach, in spirit, follows 
along similar lines as recent effort o^^'^° , but the proposed 
formulation differs importantly in the way the extended 
electrostatic interactions are treated. In particular, the 
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proposed formulation provides a unified framework to 
compute the configurational forces associated with both 
internal ionic and cell relaxations discussed in Im^ 

We begin by considering the electrostatic interactions 
that are extended in the real-space. We denote by (5(x — 
R/) a regularized Dirac distribution located at R/, and 
the nuclear charge is given by the charge distribution 
-Zi6{x - Ri). Defining p„„(x) = - X)/ Zi^\x - R/|) 
and Pn„(x') = - J2j ^J^(|x'-Rj|), the repulsive energy 
Ezz can subsequently be reformulated as 

p (x)p (x') ^^^^,_ 

|x— x'l 

where Eg^i / denotes the self energy of the nuclear charges 
and is given by 



E^ 


self 




Z/^(|x-R/|)Z/5(|x'-R/ 


-dxdx'. 

( 12 ) 

I 


We denote the electrostatic potential corresponding to 
the nuclear charge (—Z/5(|x' — R/|)) as and 

is given by 




(13) 


The self energy, thus, can be expressed as 


E. 


= -^E / ^/^'(|x-R/|)^^(x)d> 


self 


(14) 


Noting that the kernel corresponding to the extended 
electrostatic interactions in equations dH-dlsl) is the 
Green’s function of the Laplace operator, the electro¬ 
static potential and the electrostatic energy can be com¬ 
puted by taking recourse to the solution of a Poisson 
equation, or, equivalently, the following local variational 
problem: 


= J |VP^(x)| 2 dx + y'z,5(|x-R,|)P^(x)dx}, 

Vj{^)=arg J |VP^(x)|2dx + J Z/^dx - R/|)P^(x)dx} . 


(15a) 

(15b) 


In the above, iJ^(R^) denotes the Hilbert space of func¬ 
tions such that the functions and their first-order deriva¬ 
tives are square integrable on K.^. 

We next consider the electrostatic interaction energy 
corresponding to both electron and nuclear charge dis¬ 
tribution. We denote this by J{p,pnu), which is given 
by 


Jip, Pnu) — 2 


(P(x) + Pnu ))(p(xO + Pnuip^ )) 


distribution) as (/), which is given by 


<)>(x) = 


p(x') -b Pn„(x') 


dx'. 


X — X' 


(17) 


The electrostatic interaction energy of the total charge 
distribution, in terms of </>, is given by 


dbcdx.'. 


J (P, Pnu) = \ I (d(x) + PnK(x))(/)(x)dx . (18) 


(16) As before, the electrostatic interaction energy as well as 
We denote the electrostatic potential corresponding to the potential of the total charge distribution can be re- 
the total charge distribution (electron and nuclear charge formulated as the following local variational problem: 


d(p, Pnu) = - min I ^ J |V(/)(x)pdx-y(p(x)-bp„u(x))0(x)dx| 


(19a) 


(/)(x) = arp min|y- J |V(()(x)pdx - y (p(x) -b Pnu(x))(/)(x)dx| . 


(19b) 


In the above, 3^ is a suitable function space corresponding 


to the boundary conditions of the problem. In particu- 
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lar, for non-periodic problems such as isolated cluster of 
atoms y = For periodic problems, y = 

where Q denotes the unit cell and Hp^^{Q) denotes the 


space of periodic functions on Q such that the functions 
and their first-order derivatives are square integrable. 

The electrostatic interaction energy in DFT, compris¬ 
ing of Eh^ Eext and (cf. equations ©-([7])), can be 
rewritten in terms of J(p, and Egeif as 


Eh{p) + Eext{P,'R) + E^^{R) = J{p,Pnu) + X! / “ ^-^1) “ (1^ “ R-j|))p(x)dx - Eself ■ (20) 


For the sake of convenience of representation, we will 
denote by V = {F^, ..., F^} the vector containing 

the electrostatic potentials corresponding to all nuclear 
charges in the simulation domain. Using the local re- 

I 


formulation of J{p,Pnu) and Egeif (cf. equations (fT51) 
and HH)), the electrostatic interaction energy in DFT 
can now be expressed as the following local variational 
problem: 


Eh + Eext + E^z = max min Ceiicj), V, p, R) (21a) 


£ei((/', V,p, R) 


-J\V(j){x)\^dx + j (p(x)-hp„„(x))())(x)dx-l-^ J (Fpi(|x-Rj|)-F/(|x-Rj|))p(x)dx 
+ / |VF^(x)pdx-f J Z/^dx-R/|)F^(x)dx| . 

(21b) 


In the above, the minimization over F^ represents a si¬ 
multaneous minimization over all electrostatic potentials 
corresponding to I = 1,2,...,M. We note that, while 
the above reformulation of electrostatic interactions has 
been developed for pseudopotential calculations, this can 
also be extended to all-electron calculations in a straight¬ 
forward manner by using iW = (W and Zj to be the total 
nuclear charge in the above expressions. Thus, this lo¬ 
cal reformulation provides a unified framework for both 
pseudopotential as well as all-electron DFT calculations. 


We now consider the local reformulation of the ex¬ 
tended interactions in the kernel energy component of 
the WGC orbital-free kinetic energy functional (cf. ©). 
Here we adopt the recently developed local real-space re¬ 
formulation of the kernel energy^Si^i, and recall the key 
ideas and local reformulation for the sake of complete¬ 
ness. We present the local reformulation of Kq and the 
local reformulations for other kernels {Ki, Ku, K 12 ) fol¬ 
lows along similar lines. Consider the kernel energy cor¬ 
responding to Kq given by 

Tko{p) = Cf j j p“(x)i4:o(|x-x'|)p'^(x')dxdx'. 

( 22 ) 


We define potentials and given by 

^^°(x) = y is:o(|x-x'|)p“(x')dx', 

^p(x) = y ^o(|x-x'|)p'^(x')dx'. (23) 

Taking the Fourier transform of the above expressions we 
obtain 

?;0(k) =/Fo(|k|)^(k), 

;o(k) = Fo(|k|);^(k). (24) 

Following the ideas developed by Choly & Kaxiras^, Kq 
can be approximated to very good accuracy by using a 
sum of partial fractions of the following form 

where Aj, Bj, j = 1.. .m are constants, possibly com¬ 
plex, that are determined using a best fit approximation. 
Using this approximation and taking the inverse Fourier 
transform of equation (I24p . the potentials in equation 
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(l23ll reduce to 

m 

j=l 

m 

= XI W] ■ (26) 

i=i 

where uj^. (x) and (x) for j = 1.. .m are given by the 
following Helmholtz equations: 

+ B,u:% + A,B,pP = 0 . (27) 

We refer to these auxiliary potentials, a;° = 

{a;°^ ... and introduced in the 

local reformulation of the kernel energy as kernel poten¬ 
tials. Expressing the Helmholtz equations in a variational 
form, we reformulate Tkq in (l22l) as the following local 
variational problem in kernel potentials: 

Tko[p)= min, max £if„(a;°,a;;g,p), (28a) 


m P ^ 

CkM,ojI,p) = '^Cf[ / — Va;)(.(x) • Vw^.(x) 

+ (x) + (x)p“(x) + a;° . (x)p^(x) 

+ H,p(“+'5)(x)]dx|. 

(28b) 

The variational problem in equation (1281) represents a 
simultaneous saddle point problem on kernel potentials 
and for j = 1,... ,m. Following a similar pro¬ 
cedure, we construct the local variational reformulations 
for the kernel energies Tk^, and Tki 2 corresponding 
to kernels Ki, Ku and K 12 , respectively. We denote by 
/:ifi(wi,a;^,p), and £^,3 (wl^ p) 

the Lagrangians with respective kernel potentials corre¬ 
sponding to kernel energies of Ki, Ku and K 12 , respec¬ 
tively. We refer to the supplemental material for the 
numerical details of the approximations for each of the 
kernels used in the present work. 

Finally, using the local variational reformulations of 
the extended electrostatic and kernel energies, the prob¬ 
lem of computing the ground-state energy for a given po¬ 
sitions of atoms is given by the following local variational 
problem in electron-density, electrostatic potentials, and 
kernel potentials: 


£o(R) 


min max min max (Cp / p(x)®/^dx-|— 
Vp 6^ bey ey ey I J 2 

+ p) 


|V\/pWPc^x + J exc(p)p(x) dx 
min Cei((ii,V, p.R)} . 


(29) 


In the above, s denotes the index corresponding to a 
kernel, and X and y are suitable function spaces cor¬ 
responding to the boundary conditions of the problem. 
In particular, for periodic problems, y = and 

^ = {\/p\\/P S p = N}. It is convenient to 

use the substitution u(x) = -\/p(^, and enforce the in¬ 


tegral constraint in X using a Lagrange multiplier. Also, 
for the sake of notational simplicity, we will denote by uia 
and ujp the array of kernel potentials {w°, 
and {w^,respectively. Subsequently, the 
variational problem in equation (1291) can be expressed as 


Fo(R) = minmax min max £(m, Wq,, wp; R) subject to : / u'^{x)dx = N, 

u^y (p^y uj% .^y uj% / 

C{u,(t),LOa,^^p', R) = c{u) -I- (Wa, a;p,M^) -I- £c(m, A) -I- mn^^^£e;(0, V,it^,R), 

C{u) = Cp J u^°/^(x) dx-I- i J |VM(x)pdx-f J exc('a^)u^(x) dx , 

CK{u}a,U}p,U^) = Wp,M^) , 

s 

Cc{u, A) = A ■u^(x) dx — . 


( 30 ) 








B. Configurational forces 

We now turn our attention to the configurational 
forces corresponding to geometry optimization. To this 
end, we employ the approach of inner variations, where 
we evaluate the generalized forces corresponding to per¬ 
turbations of underlying space, which provides a uni¬ 
fied expression for the generalized force corresponding 
to the geometry of the simulation cell—internal atomic 
positions, as well as, the external cell domain. We con¬ 
sider infinitesimal perturbations of the underlying space 
corresponding to a generator r(x) given 
by r = such that '00 = I- We constrain the 

generator F such that it only admits rigid body deforma¬ 
tions in the compact support of the regularized nuclear 
charge distribution in order to preserve the integral 
constraint f S(x — Rj)dx = 1 . Let x denote a point 
in Q, whose image in Q' = tpeiQ) is x' = 0e(x). The 
ground-state energy on Q' is given by 

^0 ( 0 e) — I^e) (iil) 

where Ue, </>£, and are solutions of the sad¬ 
dle point variational problem given by equation dSOD 
evaluated over the function space y' = Hp^^{Q'). 
The subscript e on £ is used to denote that the 
variational problem is solved on Q' = For 

the sake of convenience, we will represent the in¬ 


tegrand of the Lagrangian £ in equation (l30l) by 
f{u, V-u, 0, V0, W/ 3 , Vw/j; Vps,Vg, R) and 

g{V~ ,\7V~ 'll), where / denotes the integrand whose 
integrals are over Q and g denotes the integrand whose 
integrals are over R^. The ground-state energy on Q' in 
terms of / and g can be expressed as 

Eoii^e) = [ /('«e(x'), Vx''Ue(x'), 00 x'), Vx' 0 e(x'),W„e(x'), 
JQ' 

Vx' e (x'), u ;/3 0x'), Vx'W/J ^ (x'); Vp0x'), (x'), 00R)) dx' 

+ E/ 5(b^'(x'),Vx-V/(x');0.(R))dx'. (32) 

j Jr3 

Transforming the above integral to domain Q, we obtain 
r 

Eo{tpe) = J /(Ue( 0 e(x)), Vx'«e( 0 e(x)). —, 0400x)), 

Vx 0 e( 0 e(x)).^,Wae( 0 £(x)), Vxa;ae( 0 e(x)). —,W^^( 0 ,(x)), 

r9x — f9x^ 

Vxu;/ 3 ,(0£(x)). —; Vps (0^(x)), (04x)), 00R)) det (—) dx 

+ E (^<^ (^)) > (0£ (x)). ^; 00R)) det (—) dx 

(33) 

We now evaluate the configurational force given by the 
Gateaux derivative of £o(0e): 


dEo{'ipe) 


de 


e =0 


f — d 5x^ 

= / (mo (x) , Vuo (x) , 00 (x) , V0O (x), Wao (x) , Vwa 0 (x), W;3o (x), (x) I Vp^ (x), (x), R) — (det (—)) 


c?x 


e =0 


+ (viV(Ix-RjI)- vq''(|x-Rj|)).(^0^ 


dx 


d0e(R,7) 


e =0 

dg 




In the above, we denote by ‘ 0 ’ the outer product be¬ 
tween two vector, by the dot product between two 
vectors and by the dot product between two tensors. 
We note that in the above expression there are no terms 
involving the explicit derivatives of / and g with respect 
to R as 5(|x' — 0e(R)|) = 5(|x — R|), which follows from 
the restriction that 0 e corresponds to rigid body defor¬ 
mations in the compact support of pnu- We further note 
that terms arising from the inner variations of £'o( 0 e) 
with respect to u^, 0e, oJae: and V~ vanish as uq 0o, 

u!ao, ^/ 3 q and V~ are the solutions of the saddle point 
variational problem corresponding to £'o(0o)- We now 


note the following identities 


d 

\ dxi 1 

dx. 

( d 

f d^k \ 

\ dxi 


de 

[dx'^j 

e =0 dx'f. 

^de 

. dxi J 

)dx'^ 

e^O 


dxj ’ 


(35) 


de 





(36) 
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Using these identities in equation (j34l) . and rearranging 
terms, we arrive at 


d£'o('0e 


de 


e =0 


= / E : Vr(x) dx + ^ / Ej : Vr(x) dx 


+ u2(x)(v(U/,-V^‘^)).(r(x)-r(Rj))dx (37) 

T 'JO 


where E and E' denote Eshelby tensors corresponding to 
/ and g, respectively. The expressions for the Eshelby 
tensors E and Ej explicitly in terms of u, </>, Wq,, Vps 
and are given by 


E = + Xu^ - ^ (U/, 


«' + E , Vu;^, I 


(38) 


(39) 


In the above, for the sake of brevity, we represented by 
fxs the integrand corresponding to Cks- We also note 
that the terms (f>Pnu and V~6{x — R/) do not appear in 

the expressions for E and Ej, respectively, as V.T = 0 
on the compact support of owing to the restriction 
that r corresponds to rigid body deformations in these 
regions. It may appear that evaluation of the second 
term in equation (13711 is not tractable as it involves an 
integral over R^. To this end, we split this integral on 
a bounded domain U containing the compact support of 
i5(x — R/), and its complement. The integral on 
can be computed as a surface integral. Thus, 

I E^ : Vr dx = / E^ : Vr dx + / e'^ : VT dx 

Jr3 Jq jR^/n 

= / Ej : Vr dx — / Ej : n 0 r ds , (40) 

Jq Jan 

where n denotes the outward normal to the surface dil. 
The last equality follows from the fact that = 0 on 

rVu. 

The configurational force in equation (1571) provides the 
generalized variational force with respect to both the in¬ 
ternal positions of atoms as well as the external cell do¬ 
main. In order to compute the force on any given atom, 
we restrict the compact support of T to only include the 
atom of interest. In order to compute the stresses asso¬ 
ciated with cell relaxation (keeping the fractional coor¬ 
dinates of atoms fixed), we restrict T to affine deforma¬ 
tions. Thus, this provides a unified expression for geom¬ 
etry optimization corresponding to both internal ionic 
relaxations as well as cell relaxation. We further note 
that, while we derived the configurational force for the 
case of pseudopotential calculations, the derived expres¬ 
sion is equally applicable for all-electron calculations by 
using = UU- 


C. Finite-element discretization 

Among numerical discretization techniques, the plane- 
wave discretization has been the most popular and widely 
used in orbital-free DFT— as it naturally lends itself 
to the evaluation of the extended interactions in electro¬ 
static energy and kernel kinetic energy functionals using 
Fourier transforms. Further, the plane wave basis offers 
systematic convergence with exponential convergence in 
the number of basis functions. However, as noted pre¬ 
viously, the plane-wave basis also suffers from notable 
drawbacks. Importantly, plane-wave discretization is re¬ 
stricted to periodic geometries and boundary conditions 
which introduces a significant limitation, especially in the 
study of defects in bulk materials^^. Further, the plane- 
wave basis has a uniform spatial resolution, and thus is 
not amenable to adaptive coarse-graining. Moreover, the 
use of plane-wave discretization involves the numerical 
evaluation of Fourier transforms whose scalability is lim¬ 
ited on parallel computing platforms. 

In order to circumvent these limitations of the plane- 
wave basis, there is an increasing focus on developing 
real-space discretization techniques for orbital-free DFT 
based on finite-difference discretization^^Ti^i and finite- 
element discretization^Si^. In particular, the finite- 
element basis^, which is a piecewise continuous poly¬ 
nomial basis, has many features of a desirable basis in 
electronic structure calculations. While being a com¬ 
plete basis, the finite-element basis naturally allows for 
the consideration of complex geometries and boundary 
conditions, is amenable to unstructured coarse-graining, 
and exhibits good scalability on massively parallel com¬ 
puting platforms. Moreover, the adaptive nature of the 
finite-element discretization also enables the considera¬ 
tion of all-electron orbital-free DFT calculations that are 
widely used in studies of warm dense matter— diiiS. Fur¬ 
ther, recent numerical studies have shown that by using a 
higher-order finite-element discretization significant com¬ 
putational savings can be realized for both orbital-free 
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DFT— and Kohn-Sham DFT calculations^!^, effectively 
overcoming the degree of freedom disadvantage of the 
finite-element basis in comparison to the plane-wave ba¬ 
sis. 

Let denote the finite-element subspace of y, where 
h represents the finite-element mesh size. The discrete 
problem of computing the ground-state energy for a given 
positions of atoms, corresponding to equation (1501) . is 
given by the constrained variational problem: 

f;o(R) = min max min max £(u/i,, cj/?. ;R) 

eyhui^g. &yh h hh ; 

3 h ^3 h 

subject to : J uf^{:>c) dx = N . (41) 

In the above, Uh, 4>h, ^ah s-nd denote the finite- 
element discretized fields corresponding to square-root 
electron-density, electrostatic potential, and kernel po¬ 
tentials, respectively. We restrict our finite-element dis¬ 
cretization such that atoms are located on the nodes 
of the finite-element mesh. In order to compute the 
finite-element discretized solution of , we represent 

S(x — Rj) as a point charge on the finite-element node 
located at Rj, and the finite-element discretization pro¬ 
vides a regularization for Previous investigations 

have suggested that such an approach provides optimal 
rates of convergence of the ground-state energy (cfi^i^ 
for a discussion). 

The finite-element basis functions also provide the gen¬ 
erator of the deformations of the underlying space in the 
isoparametric formulation, where the same finite-element 
shape functions are used to discretize both the spatial 
domain as well as the fields prescribed over the domain. 
Thus, the configurational force associated with the loca¬ 
tion of any node in the finite-element mesh can be com¬ 
puted by substituting for F, in equation (1571) . the finite- 
element shape function associated with the node. Thus, 
the configurational force on any finite-element node lo¬ 
cated at an atom location corresponds to the variational 
ionic force, which are used to drive the internal atomic 
relaxation. The forces on the finite-element nodes that 
do not correspond to an atom location represent the gen¬ 
eralized force of the energy with respect to the location of 
the finite-element nodes, and these can be used to obtain 
the optimal location of the finite-element nodes—a basis 
adaptation technique. 

We note that the local real-space variational formula¬ 
tion in section Fill A[ where the extended interactions in 
the electrostatic energy and kernel functionals are refor¬ 
mulated as local variational problems, is essential for the 
finite-element discretization of the formulation. 


IV. NUMERICAL IMPLEMENTATION 

In this section, we present the details of the numeri¬ 
cal implementation of the finite-element discretization of 
the real-space formulation of orbital-free DFT discussed 


in section imi Subsequently, we discuss the solution pro¬ 
cedure for the resulting discrete coupled equations in 
square-root electron-density, electrostatic potential and 
kernel potentials. 

A. Finite-element basis 

A finite-element discretization using linear tetrahedral 
finite-elements has been the most widely used discretiza¬ 
tion technique for a wide range of partial differential 
equations. Linear tetrahedral elements are well suited 
for problems involving complex geometries and moder¬ 
ate levels of accuracy. However in electronic structure 
calculations, where the desired accuracy is commensu¬ 
rate with chemical accuracy, linear finite elements are 
computationally inefficient requiring of the order of hun¬ 
dred thousand basis functions per atom to achieve chem¬ 
ical accuracy. A recent studj*2i has demonstrated the 
significant computational savings—of the order of 1000- 
fold compared to linear finite-elements—that can be real¬ 
ized by using higher-order finite-element discretizations. 
Thus, in the present work we use higher-order hexahedral 
finite elements, where the basis functions are constructed 
as a tensor product of basis functions in one-dimension^. 

B. Solution procedure 

The discrete variational problem in equation (1411) in¬ 
volves the computation of the following fields—square- 
root electron-density, electrostatic potential and kernel 
potentials. Two solution procedures, suggested in prior 
efforts^, for solving this discrete variational problem in¬ 
clude: (i) a simultaneous solution of all the discrete fields 
in the problem; (ii) a nested solution procedure, where 
for every trial square-root electron-density the discrete 
electrostatic and kernel potential fields are computed. 
Given the non-linear nature of the problem, the simul¬ 
taneous approach is very sensitive to the starting guess 
and often suffers from lack of robust convergence, espe¬ 
cially for large-scale problems. The nested solution ap¬ 
proach, on the other hand, while constituting a robust 
solution procedure, is computationally inefficient due to 
the huge computational costs incurred in computing the 
kernel potentials which involves the solution of a series 
of Helmholtz equations (cf. equation (157)) b Thus, in the 
present work, we will recast the local variational prob¬ 
lem in equation (1411) as the following fixed point iteration 
problem: 

{uh,M = arg min arg uiaxC{uh, 4>h,3Jag,,3ji3^-,'R) 

rph 

subject to : / 'u^(x) dx = N. (42a) 

= arg min arg max£(«?;, , W/ 3 ;.; R) • 

(42 b) 

We solve this fixed point iteration problem using a mix¬ 
ing scheme, and, in particular, we employ the Anderson 




11 


mixing scheme^ with full history in this work. Our nu¬ 
merical investigations suggest that the fixed point itera¬ 
tion converges, typically, in less than ten self-consistent 
iterations even for large-scale problems, thus, provid¬ 
ing a numerically efficient and robust solution procedure 
for the solution of the local variational orbital-free DFT 
problem. We note that this idea of fixed point iteration 
has independently and simultaneously been investigated 
by another group in the context of finite difference dis¬ 
cretization^, and have resulted in similar findings. 

In the fixed point iteration problem, we employ a si¬ 
multaneous solution procedure to solve the non-linear 
saddle point variational problem in Uh and (j)^ (equa¬ 
tion (|42all l. We employ an inexact Newton solver pro¬ 
vided by the PETSc package^ with held split precondi¬ 
tioning and generalized-minimal residual method (GM- 
RES)^ as the linear solver. The discrete Helmholtz equa¬ 
tions in equation (j42bp are solved by employing block 
Jacobi preconditioning and using GMRES as the linear 
solver. An efhcient and scalable parallel implementation 
of the solution procedure has been developed to take ad¬ 
vantage of the parallel computing resources for conduct¬ 
ing the large-scale simulations reported in this work. 


V. RESULTS AND DISCUSSION 

In this section, we discuss the numerical studies on 
Al, Mg and Al-Mg intermetallics to investigate the 
accuracy and transferability of the real-space formu¬ 
lation of orbital-free DFT (RS-OFDFT) proposed in 
section nni Wherever applicable, we benchmark the 
real-space orbital-free DFT calculations with plane-wave 
based orbital-free DFT calculations conducted using 
PROFESS^l, and compare with Kohn-Sham DFT (KS- 
DFT) calculations conducted using the plane-wave based 
ABINIT code^i^. Further, we demonstrate the useful¬ 
ness of the proposed real-space formulation in studying 
the electronic structure of isolated defects. 


A. General calculation details 

In all the real-space orbital-free DFT calculations re¬ 
ported in this section, we use the local reformulation of 
the density-dependent WGGii kinetic energy functional 
proposed in section IIII Al the local density approxima¬ 
tion (LDA)^ for the exchange-correlation energy, and 
bulk derived local pseudopotentials (BLPS)^ for Al and 
Mg. Gell stresses and ionic forces are calculated using the 
unified variational formulation of configurational forces 
developed in section llilBl In the second order Taylor ex¬ 
pansion of the density-dependent WGG functional about 
the bulk electron density (cf. Section |n|, we only retain 
the Ki 2 term for the computation of bulk properties as 
the contributions from K 12 dominate those of Kn for 
bulk materials systems. However, in the calculations in¬ 
volving mono-vacancies, where significant spatial pertur- 



FIG. I: Convergence of the finite-element approximation in 
the energy of a fee Al unit cell with lattice constant a = 7.2 

Bohr. 



FIG. 2: Convergence of the finite-element approximation in 
the hydrostatic stress of a fee Al unit cell with lattice 
constant a = 7.2 Bohr. 


bations in the electronic structure are present, we use the 
full second order Taylor expansion of the density depen¬ 
dent WGG functional. We recall from section PlI Al that 
in order to obtain a local real-space reformulation of the 
extended interactions in the kinetic energy functionals, 
the kernels (ATo, Ki, Kn, K 12 ) are approximated using 
a sum of m partial fractions where the coefficients of the 
partial fractions are computed using a best fit approxima¬ 
tion (cf. equation (l25ll b These best fit approximations 
for m = 4,5,6 that are employed in the present work are 
given in the supplemental material. It has been shown in 
recent studies that m = 4 suffices for Al^i^. However, 
we find that m = 6 is required to obtain the desired ac¬ 
curacy in the bulk properties of Mg, and Table HIl shows 
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the comparison between the kernel approximation with 
m = 6 and plane-wave based orbital-free DFT calcula¬ 
tions conducted using PROFESS^ for Mg. Thus, we use 
the best fit approximation of the kernels with m = 4 for 
Al, and employ the approximation with m = 6 for Mg 
and Al-Mg intermetallics. Henceforth, we will refer by 
RS-OFDFT-FE the real-space orbital-DFT calculations 
conducted by employing the local formulation and finite 
element discretization proposed in section imi 

The KS-DFT calculations used to assess the accuracy 
and transferability of the proposed real-space orbital- 
free DFT formalism are performed using the EDA ex¬ 
change correlation functional^. The KS-DFT calcu¬ 
lations are conducted using both local BLPS as well 
as the non-local Troullier-Martins pseudopotential (TM- 
NLPS)^ in order to assess the accuracy and transfer- 
ability of both the model kinetic energy functionals in 
orbital-free DFT as well as the local pseudopotentials to 
which the orbital-free DFT formalism is restricted to. 
The TM-NLPS for Al and Mg are generated using the 
fhi98PP code^. Within the fhi98PP code, we use the fol¬ 
lowing inputs: 3d angular momentum channel as the local 
pseudopotential component for both Al and Mg, default 
core cutoff radii for the 3s, 3p, and 3d angular momen¬ 
tum channels, which are {1.790, 1.974, 2.124} Bohr and 
(2.087, 2.476, 2.476} Bohr for Al and Mg respectively, 
and the LDA^l exchange-correlation. For brevity, hence¬ 
forth, we refer to the KS-DFT calculations with BLPS 
and TM-NLPS as KS-BLPS and KS-NLPS, respectively. 

In all the RS-OFDFT-FE calculations reported in this 
work, the finite-element discretization, order of the finite- 
elements, numerical quadrature rules and stopping tol¬ 
erances are chosen such that we obtain 1 meV/ atom 
accuracy in energies, 1 x 10~^ Hartree Bohr“^ accuracy 
in cell stresses and 1 x 10“® Hartree Bohr”^ accuracy 
in ionic forces. Similar accuracies in energies, stresses 
and ionic forces are achieved for KS-DET calculations 
by choosing the appropriate k-point mesh, plane-wave 
energy cutoff, and stopping tolerances within ABINIT’s 
framework. All calculations involving geometry opti¬ 
mization are conducted until cell stresses and ionic forces 
are below threshold values of 5 x 10“^ Hartree Bohr“^ 
and 5 x 10“^Hartree Bohr“^, respectively. 


B. Convergence of finite-element discretization 

We now study the convergence of energy and stresses 
with respect to the finite-element discretization of the 
proposed real-space orbital-free DFT formulation. In a 
prior study on the computational efficiency afforded by 
higher-order finite-element discretization in orbital-free 
DFT—, it was shown that second and third-order finite- 
elements offer an optimal choice between accuracy and 
computational efficiency. Thus, in the present study, 
we limit our convergence studies to HEX27 and HEX64 
finite-elements, which correspond to second- and third- 
order finite-elements. As a benchmark system, we con¬ 
sider a stressed fee Al unit cell with a lattice constant 


a = 7.2 Bohr. We first construct a coarse finite-element 
mesh and subsequently perform a uniform subdivision to 
obtain a sequence of increasingly refined meshes. We de¬ 
note by h the measure of the size of the finite-element. 
For these sequence of meshes, we hold the cell geometry 
fixed and compute the discrete ground-state energy, Eh, 
and hydrostatic stress, Oh- The extrapolation procedure 
proposed in Motamarri et. al^i allows us to estimate the 
ground-state energy and hydrostatic stress in the limit 
as h —>■ 0, denoted by Eq and ctq. To this end, the en¬ 
ergy and hydrostatic stress computed from the sequence 
of meshes using HEX64 finite-elements are fitted to ex¬ 
pressions of the form 

\cro - ah\ = Ca , (43) 

to determine Eq, qe, (Tq, Sz qa-. In the above expression, 
Nei denotes the number of elements in a finite-element 
mesh. We subsequently use Eq and ag as the exact val¬ 
ues of the ground-state energy and hydrostatic stress, 
respectively, for the benchmark system. Figures [T] and [2] 

show the relative errors in energy and hydrostatic stress 
1 

plotted against which represents a measure of 

h. We note that the slopes of these curves provide the 
rates of convergence of the finite-element approximation 
for energy and stresses. These results show that we ob¬ 
tain close to optimal rates of convergence in energy of 
where k is polynomial interpolation order {k = 2 
for HEX27 and fc = 3 for HEX64). Further, we obtain 
close to convergence in the stresses, which rep¬ 

resents optimal convergence for stresses. The results also 
suggest that higher accuracies in energy and stress are 
obtained with HEX64 in comparison to HEX27. Thus, 
we will employ HEX64 finite-elements for the remainder 
of our study. 


TABLE I: The energy difference in eV between a stable 
phase and the most stable phase for Al and Mg computed 
using RS-OFDFT-FE and KS-DFT with TM-NLPS. 


Al 

fee 

hep 

bcc 

SC 

dia 

RS-OFDFT-FE 

qa 

0.016 

0.075 

0.339 

0.843 

KS-NLPS 

0 

0.038 

0.106 

0.400 

0.819 

Mg 

hep 

fee 

bcc 

SC 

dia 

RS-OFDFT-EE 

0 

0.003 

0.019 

0.343 

0.847 

KS-NLPS 

0 

0.014 

0.030 

0.400 

0.822 


^ The zero in the first column is to indicate that these numbers 
are the reference against which energies of other phases are 
determined. 
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TABLE II: Bulk properties of A1 and Mg: Equilibrium ground-state energy per atom (Emin in eV), volume per atom (Vb in 
A®) and bulk modulus (Bo in GPa) computed using RS-OFDFT-FE, PROFESS, and KS-DFT with BLPS and TM-NLPS. 


AP 

RS-OFDFT-FE 

PROFESS 

KS-BLPS 

KS-NLPS 

E'min 

-57.935 

-57.936 

-57.954 

-57.207 

Vo 

15.68 

15.68 

15.62 

15.55 

Bo 

81.7 

81.5 

84.1 

83.6 

Mg^ 

RS-OFDFT-FE 

PROFESS 

KS-BLPS 

KS-NLPS 

-S^min 

-24.647 

-24.647 

-24.678 

-24.514 

K) 

21.40 

21.43 

21.18 

21.26 

Bo 

36.8 

36.6 

38.5 

38.6 


^ Cell-relaxed lattice constant tor fee A1 using RS-OFDFT-FE is 
ao = 7.51 Bohr. 

Cell-relaxed lattice constants for hep Mg using RS-OFDFT-FE 
are ag = 5.89 Bohr, eg = 9.62 Bohr. 


C. Bulk properties of Al, Mg and Al-Mg 
intermetallics 

We now study the accuracy and transferability of the 
proposed real-space formulation of orbital-free DFT for 
bulk properties of Al, Mg and Al-Mg intermetallics. To 
this end, we begin with the phase stability study of Al 
and Mg, where we compute the difference in the ground- 
state energy of a stable phase and the ground-state en¬ 
ergy of the most stable phase. The results for Al and Mg 
are shown in Table HI and are compared against those ob¬ 
tained with KS-DFT employing TM-NLPS. We note that 
RS-OFDFT-FE correctly predicts the most stable phases 
of Al and Mg being fee and hep, respectively. Further, the 
stability ordering of the various phases computed using 
RS-OFDFT-FE is consistent with KS-DFT TM-NLPS 
calculations. Moreover, the energy differences between 
the various stable phases and the most stable phase com¬ 
puted using RS-OFDFT-FE are in close agreement with 
KS-DFT calculations. 

We next consider bulk properties of Al, Mg and Al-Mg 
intermetallics. To this end, for each system, we first op¬ 
timize cell geometry and ionic positions to determine the 
equilibrium cell structure, equilibrium volume (Vq) a-nd 
ground-state energy (Emin). We subsequently compute 
the bulk modulus given by^ 


B = V 


d^E 


dV^ 


V=Vo 


(44) 


where E denotes the ground-state energy of a unit-cell 
with volume V. To compute the bulk modulus, we vary 
the cell volume by applying a volumetric deformation to 
the relaxed (equilibrium) unit-cell, which transforms the 
equilibrium cell vectors {ci , C 2 , 03 } to {c^ ,C 2 ,€ 3 } and 
are given by 


cb = cy (1 + v). (45) 

While keeping the cell structure fixed, we calculate the 


ground-state energy for each 77 between — 0.01 to 0.01 in 
steps of 0.002 and fit a cubic polynomial to the E — V 
data. We subsequently compute the bulk modulus, using 
equation (l44l) . at the equilibrium volume, Vq. The com¬ 
puted bulk properties—ground-state energy, equilibrium 
volume and bulk modulus at equilibrium—for Al and Mg 
are given in Table im and those of Al-Mg intermetallics 
(Al 3 Mg, Mgi 3 Ali 4 , Mg;i 7 Ali 2 , and Mg23Al3o) are given 
in Table uni These results suggest that the bulk proper¬ 
ties of Al, Mg and Al-Mg intermetallics computed using 
RS-OFDFT-FE are in good agreement with PROFESS 
and KS-DFT calculations. 

Finally, we consider the formation energies of Al-Mg 
intermetallics. In addition to the Al-Mg intermetallics 
for which we computed the bulk properties, we also com¬ 
pute the formation energy of the /3' alloy. The /3' alloy 
has a disorder in 10 out of 879 sites with each site having 
0.5 chance of being occupied by either Al or Mg^. In 
our simulations, we consider the two limits where all 10 
sites are occupied by either Al or Mg and refer to these as 
/3'(A1) and /3'(Mg), respectively. For these two systems, 
we do not provide KS-DFT results as they are computa¬ 
tionally prohibitive. The formation energies for the range 
of Al-Mg intermetallics are reported in Table ITVl Our re¬ 
sults suggest that the formation energies predicted by 
RS-OFDFT-FE are in good agreement with PROFESS 
calculations, and in close agreement with KS-DFT cal¬ 
culations. 


D. Configurational forces and atomic displacements 

As a next step in our study of the accuracy and trans¬ 
ferability of RS-OFDFT-FE, we compute the configura¬ 
tional forces on atoms that are perturbed from their equi¬ 
librium positions and compare these with Kohn-Sham 
DFT calculations. We investigate the accuracy of the 
forces in both fee Al and hep Mg. We begin by consider¬ 
ing the relaxed Al fee unit cell, and the relaxed Mg hep 
unit cell. In the relaxed Al fee unit cell, we perturb the 
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TABLE III: Bulk properties of Al-Mg intermetallics: Equilibrium ground-state energy per primitive cell (Amin in eV), 
volume of primitive cell (Vo in A^), and bulk modulus (Bo in GPa) computed using RS-OFDFT-FE, PROFESS, and KS-DFT 

with BLPS and TM-NLPS. 


AlaMg 

RS-OFDFT-FE 

PROFESS 

KS-BLPS 

KS-NLPS 

E • 

-^min 

-198.492 

-198.496 

-198.575 

-I96.I62 

^0 

67.23 

67.31 

67.13 

66.52 

Bo 

69.2 

67.0 

67.6 

71.0 

Mgi3Ali4 

RS-OFDFT-FE 

PROFESS 

KS-BLPS 

KS-NLPS 

E • 

-^min 

-II30.083 

-II30.I00 

-II30.972 

-III7.936 


494.77 

494.73 

498.19 

492.73 

Bo 

53.1 

52.1 

54.7 

54.8 

Mg47Ali2 

RS-OFDFT-FE 

PROFESS 

KS-BLPS 

KS-NLPS 

-S^min 

-III4.446 

-III4.526 

-III6.I85 

-II04.0I2 

Ko 

545.32 

544.85 

543.67 

544.21 

Bo 

51.1 

52.3 

55.2 

54.4 

MgssAlao 

RS-OFDFT-FE 

PROFESS 

KS-BLPS 

KS-NLPS 

Emin 

-2306.785 

-2306.762 

-2307.989 

-2281.082 

^0 

953.87 

952.55 

963.72 

957.46 

Bo 

64.2 

60.9 

60.5 

60.5 


TABLE IV: Formation energy per atom (eV/atom) of Al-Mg intermetallics calculated using RS-OFDFT-FE, PROFESS, 

and KS-DFT with TM-NLPS. 


Method 

AlaMg 

Mg43Ali4 

Mgl7Ali2 

Atg23Al30 

^'(Al) 

/3'(Mg) 

RS-OFDFT-FE 

- 0.010 

0.053 

-0.008 

-0.035 

-0.026 

- 0.020 

PROFESS 

-O.OII 

0.052 

-O.OII 

-0.034 

-0.029 

-0.023 

KS-NLPS 

-0.007 

0.061 

-0.027 

-0.019 

- 

- 


face-centered atom with fractional coordinates 0 , ^ 

by 0.1 Bohr in the [0 1 0] direction. In the relaxed Mg hep 
unit cell, we perturb the atom with fractional coordinates 
|, i, i by 0.1 Bohr in the [2 I 3 0] direction (directions 
in hep Mg are represented using Miller-Bravais indices). 
The conhgurational forces on the perturbed atoms are 
computed using RS-OFDFT-FE, and compared against 
KS-DFT calculations. The computed restoring forces, 
along [0 1 0] for the A1 system and along [2 13 0] for 
the Mg system, are reported in Table El We note that 
the computed restoring forces from RS-OFDFT-FE are 
in good agreement with PROFESS and KS-DFT calcu¬ 
lations. 

As a more stringent test of accuracy and transferabil¬ 
ity, we consider the atomic relaxations around a mono¬ 
vacancy in fee A1 and hep Mg. In the case of mono¬ 
vacancy in Al, we consider a supercell containing 3x3x3 
fee Al unit cells and remove an atom to create a mono¬ 
vacancy. We calculate the forces on the neighboring 
atoms of the mono-vacancy, and their relaxation displace¬ 
ments upon ionic relaxation using both RS-OFDFT-FE 
and KS-DFT calculations. Periodic boundary conditions 


are employed in these calculations. Table Ell reports the 
computed force and relaxation displacement in Al on the 
nearest neighboring atom, which experiences the largest 
ionic force and relaxation. In the case of a mono-vacancy 
in Mg, we consider a supercell containing 3x3x2 hep 
unit cells, and Table lVlIl reports the ionic force and relax¬ 
ation displacement on the neighboring atom that has the 
largest force in the presence of the vacancy. As is evident 
from the results, the ionic forces and relaxed displace¬ 
ments for a mono-vacancy in Al and Mg computed using 
RS-OFDFT-FE are in good agreement with PROFESS, 
and in close agreement with KS-DFT calculations. These 
results suggest that the proposed real-space orbital-free 
DFT formulation provides a good approximation to KS- 
DFT for Al-Mg materials systems. 


E. Cell-size studies on a mono-vacancy in Al 

Prior Fourier-space calculations using OF-DFT and 
WGC Functional^^, and KS-DFT calculations^ have 
suggested that cell-sizes containing ^ 256 lattice sites 
are sufficient to obtain a well-converged (to within 3 
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TABLE V: Restoring force (eV/Bohr) on the perturbed 
atom in fee A1 and hep Mg unit cells computed using 
RS-OFDFT-FE, PROFESS, and KS-DFT calculations. 



RS-OFDFT-FE 

PROFESS 

KS-BLPS 

KS-NLPS 

Al 

0.148 

0.137 

0.134 

0.126 

Mg 

0.019 

0.019 

0.018 

0.019 


TABLE VI: Ionic forces (eV/Bohr) and relaxation 
displacement (Bohr) on the nearest neighboring atom to a 
mono-vacancy in a periodic 3x3x3 fee A1 supercell, 
calculated using RS-OFDFT-FE, PROFESS, and KS-DFT. 
/ and d denote the magnitudes of ionic force and relaxation 
displacement. Zf and Zd denote the angles (in degrees) of 
the force and displacement vectors with respect to the 
KS-NLPS force and displacement vectors. 



RS-OFDFT-FE 

PROFESS KS-BLPS 

KS-NLPS 

/ 

0.141 

0.146 0.130 

0.119 

d 

9.90x10-2 

9.75x10-2 9.47x10-2 

8.90x10-2 

Zf 

0.00 

0.00 0.00 

0.00 

Zd 

0.15 

0.00 0.00 

0.00 


meV) mono-vacancy formation energy in fee Al. These 
Fourier-space calculations, which employ periodic bound¬ 
ary conditions, compute the properties of a periodic ar¬ 
ray of vacancies. On the other hand, real-space calcu¬ 
lations on isolated mono-vacancies in bulk, computed 
using the recently developed coarse-graining techniques 
for orbital-free DFT— suggest that cell-size effects in 
mono-vacancy calculations are present up to cell-sizes of 
^ 10^ atoms. Although both approaches give similar 
converged vacancy formation energies, this discrepancy 
in the cell-size effects has thus far remained an open ques¬ 
tion. 

In order to understand the source of this discrepancy, 
we conduct a cell-size study of the mono-vacancy forma¬ 
tion energy in Al using RS-OFDFT-FE with two types 
of boundary conditions: (i) periodic boundary conditions 
on electronic fields; (ii) Dirichlet boundary conditions 
on electronic fields with values corresponding to that of 


TABLE VII: Ionic forces (eV/Bohr) and relaxation 
displacement (Bohr) on the nearest neighboring atom to a 
mono-vacancy in a periodic 3x3x2 hep Mg supercell, 
calculated using RS-OFDFT-FE, PROFESS, and KS-DFT. 



RS-OFDFT-FE 

PROFESS KS-BLPS 

KS-NLPS 

/ 

0.059 

0.060 0.053 

0.046 

d 

8.26x10-2 

8.64x10-2 7.00x10-2 

5.83x10-2 

Zf 

5.11 

4.73 2.75 

0.0 

Zd 

5.66 

5.27 3.58 

0.0 


a perfect crystal. These Dirichlet boundary conditions, 
which we refer to as bulk Dirichlet boundary conditions, 
correspond to the scenario where perturbations in the 
electronic structure due to the mono-vacancy vanish on 
the boundary of the computational domain, and the elec¬ 
tronic structure beyond the computational domain cor¬ 
responds to that of the bulk. We note that periodic 
boundary conditions mimic the widely used Fourier-space 
calculations on point defects, whereas the bulk Dirich¬ 
let boundary conditions correspond to simulating an iso¬ 
lated point defect embedded in bulk. We note that the 
local real-space formulation of orbital-free DFT and the 
finite-element basis are key to being able to consider these 
boundary conditions. 


We compute the vacancy formation at constant volume 

asi2^ 

/ N -1 \ N -I 

E,f = Ei^^N-1,1, ^^n]-^^EiN,o,n), 

(46) 

where E {N, 0, D) denotes the energy of perfect crys¬ 
tal containing N atoms occupying a volume D, and 
E{N — 1,1, denotes energy of a computational 

cell containing A^ — 1 atoms and one vacancy occupying 
a volume For both periodic boundary conditions 

and bulk Dirichlet boundary conditions, the lattice site 
where the vacancy is created is chosen to be the farthest 
site from the domain boundary. As we are primarily in¬ 
terested in the cell-size effects of the electronic structure, 
we do not consider ionic relaxations in this part of our 
study. Table IVIIII shows the unrelaxed mono-vacancy 
formation energies for different cell sizes computed using 
RS-OFDFT-FE using both periodic boundary conditions 
and bulk Dirichlet boundary conditions. We note that 
the mono-vacancy formation energies using both sets of 
boundary conditions converge to the same value, and this 
is also in good agreement with PROFESS and KS-DFT 
calculations (cf. Table HXl) . However, it is interesting to 
note that the mono-vacancy formation energies with pe¬ 
riodic boundary conditions are well converged (to within 
10 meV) by 3 x 3 x 3 cell-size (108 atoms), whereas we 
required a 6 x 6 x 6 cell-size (864 atoms) to achieve a 
converged formation energy with bulk Dirichlet bound¬ 
ary conditions. 


In order to understand this boundary condition depen¬ 
dence of the cell-size effects, we compute the perturba¬ 
tions in the electronic fields due to the presence of the 
mono-vacancy by subtracting from the electronic fields 
corresponding to the mono-vacancy the electronic fields 
of a perfect crystal. To this end, we define the normal¬ 
ized perturbations in the electronic fields computed on 
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TABLE VIII: Unrelaxed mono-vacancy formation energies 
for A1 computed using RS-OFDFT-FE with periodic 
boundary conditions (Ajj- in eV) and bulk Dirichlet 
boundary conditions in eV). 


Cell size 

N 

j^bD 

Kf 

2 x 2 x 2 

32 

-0.390 

0.955 

3x3x3 

108 

0.864 

0.915 

4x4x4 

256 

0.971 

0.908 

5x5x5 

500 

0.944 

- 

6 x 6 x 6 

864 

0.918 

- 

7x7x7 

1372 

0.914 

- 


TABLE IX: Unrelaxed mono-vacancy formation energies 
[E^f in eV) for Al computed using PROFESS^, and 
KS-DFT on a 3 X 3 X 3 computational cell. 


Eyf 

PROFESS 

0.903 

KS-DFT-BLPS 

0.815 

KS-DFT-NLPS 

0.811 


the finite-element mesh to be 

K = (Uh - O /Vav «) , 
/Vav ^ 


K,h = " Z! /Vav 


,h ’ 




i=i 


vi=i 




P 
0j,h 


/Vav 


P 
0j A 




i=i 




(47) 


In the above, {uh,4>h,uJaj,h,uJi3j,h} and 
h} denote the electronic fields in 
the computational domain with the vacancy and those 
without the vacancy (perfect crystal), respectively. 
Vav(.) denotes the volume average of an electronic field 
over the computational cell. As a representative metric, 
in the definition of ^ and we only consider the 
kernel potentials corresponding to Kq. Figures [ 3 ] and 0 ] 
shows the normalized corrector fields for the mono¬ 
vacancy, computed using periodic boundary conditions, 
along the face-diagonal of the periodic boundary. It is 
interesting to note from these results that the pertur¬ 
bations in the electronic structure due to the vacancy 
are significant up to 6 x 6 x 6 computational cells. 
Thus, although the vacancy formation energy appears 
converged by 3 x 3 x 3 computational cell while using 
periodic boundary conditions, the electronic fields are 
not converged till a cell-size of 6 x 6 x 6 computational 


cell. On the other hand, the cell-size convergence in 
mono-vacancy formation energy suggested by the bulk 
Dirichlet boundary conditions is inline with the conver¬ 
gence of electronic fields. These results unambiguously 
demonstrate that the cell-size effects in the electronic 
structure of defects are larger than those suggested by 
a cell-size study of defect formation energies employing 
periodic boundary conditions. Using bulk Dirichlet 
boundary conditions for the cell-size study of defect 
formation energies provides a more accurate estimate of 
the cell-size effects in the electronic structure of defects, 
and the extent of electronic structure perturbations due 
to a defect. Further, while periodic boundary conditions 
are limited to the study of point defects, bulk Dirichlet 
boundary conditions can be used to also study defects 
like isolated dislocations^^, whose geometry does not 
admit periodic boundary conditions. 


VI. SUMMARY 

We have developed a local real-space formulation of 
orbital-free DFT with WGC kinetic energy functionals by 
reformulating the extended interactions in electrostatic 
and kinetic energy functionals as local variational prob¬ 
lems in auxiliary potentials. The proposed real-space 
formulation readily extends to all-electron orbital-free 
DFT calculations that are commonly employed in warm 
dense matter calculations. Building on the proposed 
real-space formulation we have developed a unified vari¬ 
ational framework for computing configurational forces 
associated with both ionic and cell relaxations. Fur¬ 
ther, we also proposed a numerically efficient approach 
for the solution of ground-state orbital-free DFT prob¬ 
lem, by recasting the local saddle point problem in the 
electronic fields—electron density and auxiliary potential 
fields—as a fixed point iteration problem and employing 
a self-consistent iteration procedure. We have employed 
a finite-element basis for the numerical discretization of 
the proposed real-space formulation of orbital-free DFT. 
Our numerical convergence studies indicate that we ob¬ 
tain close to optimal rates of convergence in both ground- 
state energy and configurational forces with respect to 
the finite-element discretization. 

We subsequently investigated the accuracy and trans¬ 
ferability of the proposed real-space formulation of 
orbital-free DFT for Al-Mg materials system. To this 
end, we conducted a wide range of studies on Al, Mg 
and Al-Mg intermetallics, including computation of bulk 
properties for these systems, formation energies of Al-Mg 
intermetallics, and ionic forces in bulk and in the pres¬ 
ence of point defects. Our studies indicate that orbital- 
free DFT and the proposed real-space formulation is in 
good agreement with Kohn-Sham DFT calculations using 
both local pseudopotentials as well as non-local pseud- 
potentials, thus providing an alternate linear-scaling ap¬ 
proach for electronic structure studies in Al-Mg materi¬ 
als system. We finally investigated the cell-size effects 
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d d 


FIG. 3: Normalized corrector fields for a mono-vacancy, computed with periodic boundary conditions, along the face 
diagonal on the computational domain boundary. The abscissa d represents a normalized coordinate along the face diagonal. 
Results for computational cell sizes from 2x2x2 to 4x4x4 are shown. 


in the electronic structure of a mono-vacancy in Al, and 
demonstrated that the cell-size convergence in the va¬ 
cancy formation energy computed by employing periodic 
boundary conditions is not commensurate with the con¬ 
vergence of the electronic fields. On the other hand, the 
true cell-size effects in the electronic structure are re¬ 
vealed by employing the bulk Dirichlet boundary con¬ 
ditions, where the perturbations in the electronic fields 
due to the defect vanish on the boundary of the com¬ 
putational domain. Our studies indicate that the true 
cell-size effects are much larger than those suggested by 
periodic calculations even for simple defects like point 
defects. We note that the proposed real-space formula¬ 
tion and the finite-element basis are crucial to employing 
the bulk Dirichlet boundary conditions that are otherwise 
inaccessible using Fourier based formulations. The pro¬ 
posed formulation, besides being amenable to complex 
geometries, boundary conditions, and providing excellent 
scalability on parallel computing platforms, also enables 
coarse-graining techniques like the quasi-continuum re- 
ductior4^i^ to conduct large-scale electronic structure 


calculations on the energetics of extended defects in Al- 
Mg materials system, and is an important direction for 
future studies. 
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FIG. 4: Normalized corrector fields for a mono-vacancy, computed with periodic boundary conditions, along the face 
diagonal on the computational domain boundary, for cell sizes ranging from 5x5x5 to 7x7x7. 
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